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I  SCHOTTKY  BARRIERS  WORK 


The  Schottky  barriers  work  concentrated  this  year  on  transport  in  the  depletion  region. 
Transport  includes  both  the  depletion  region  and  the  quasineutral  region,  and  phenomena  caus¬ 
ing  departures  from  ideal  electrical  behavior  occurs  in  both.  In  last  year’s  report,  however,  it 
was  shown  that  the  depletion  region  can  and  should  be  treated  separately  from  the  quasineutral 
region;  the  two  can  be  coupled  together  by  effective  boundary  conditions  at  the  interface. 
What  is  required  in  the  depiction  region  are  the  scattering  probabilities  (or  equivalently  the 
transmission  and  reflection  coefficients)  for  an  incident  electron  entering  the  depletion  region. 
Last  year’s  work  centered  on  numerical  solutions  of  the  Boltzmann  transport  equations  in  the 
quasineutral  region  for  a  highly  non-Boltzmann,  spatially  varying  electron  distribution  subject 
to  boundary  conditions  at  the  interface  connecting  the  quasineutral  and  depletion  regions.  This 
year’s  work  has  focused  on  depiction  region,  to  obtain  the  scattering  rates  of  electrons  travers¬ 
ing  it.  At  this  time  last  year,  it  was  believed  that  scattering  from  a  collection  of  random 
impurities  (in  contradistinction  to  a  jellium  of  charge  which  gives  rise  to  the  usual  one¬ 
dimensional  quadratic  potential  barrier)  would  be  very  large  (Boudville  and  McGill,  1985)  and 
would  require  a  fully  three-dimensional  solution  the  time-dependent  Schrodinger  equation  for 
an  (effective  mass)  electron  moving  through  the  depletion  region.  A  careful  study  of  the  effect 
of  a  distribution  of  random  impurities  is  discussed  in  the  following  pages.  The  barrier  result¬ 
ing  from  a  distribution  of  random  impurities  was  found — somewhat  surprisingly — to  differ 
remarkably  little  from  that  of  a  jellium.  The  difference  between  these  potentials  acts  as  a  weak 
perturbation  to  electrons  traversing  the  barrier. 

This  result  greatly  simplifies  the  calculation  of  scattering  in  the  depletion  region,  because 
now  coulomb  scattering  can  possibly  be  neglected,  or  be  treated  in  perturbation  theory  on  the 
same  footing  as  the  other  scattering  mechanisms  (e.g.,  intervalley  scattering  and  optical  phonon 
scattering).  Electron  scattering  is  dominated  by  a  smooth  one-dimensional  potential,  which 
includes  the  parabolic  potential  of  jellium,  the  1/z  potential  from  image-force  lowering,  and 
possibly  another  short-range  potential  from  metal-induced  gap  states  or  intcrfacial  states. 
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We  can  exploit  the  smallness  of  the  perturbation  in  the  following  way.  The  scattering  rate 
w  is  given  by  Fermi’s  golden  rule  in  terms  of  the  T  matrix  (Schiff,  1968): 

w  =  -y  p(P)  l(piTla)l 2 

where  p(p)  is  the  density  of  outgoing  states  and  the  T  matrix 
I  (J3 1 T I  ct)  I  =  jup‘(f)V(f)^d3r 

represents  the  transmission  amplitude  coupling  a  free  electron  state  up  to  a  scattered  outgoing 
state  §a.  (£,  is  a  solution  to  the  full  Hamiltonian,  including  the  potential  V.) 

By  expressing  the  scattering  potential  as  a  sum  of  two  terms  V!  and  Vp,  where  V,  is  the 
large  one-dimensional  potential  and  Vp  is  the  small  perturbation,  it  turns  out  to  be  possible  to 
express  the  T  matrix  as  a  sum  of  a  part  that  arises  from  V,  alone  and  a  correction  term  (Schiff, 
1968).  The  T  matrix  is  given  by 

l(piTla)l  =  l(plT,la)l  + /^(OVp^^r 

=  l(J3IT,la)l  +  J^pOOVp^a^r 

In  the  approximate  expression,  the  second  term  contains  a  matrix  element  of  Vp  between  states 
distorted  by  the  Vj.  This  approximation  consists  in  replacing  in  the  second  term  the  true  scat¬ 
tered  state  by  one  distorted  only  by  Vp  Because  Vp  is  small,  the  approximation  should  be 
quite  good,  essentially  exact  for  our  purposes.  Because  Vt  is  only  one -dimensional,  it  is  tract¬ 
able  to  obtain  the  T  matrix  exactly  for  Vt,  the  resulting  wave  functions  ^  can  be  used  to  obtain 
the  second  terms  in  the  above  expression  and  the  total  scattering  rate  can  be  evaluated. 

To  this  end,  we  have  developed  a  code  that  solves  numerically  the  wave  functions  and  T 
matrix  for  effective  mass  electrons  in  an  arbitrary  one-dimensional  potential.  There  are  some 
difficulties  associated  with  numerical  solutions,  because  the  wave  functions  have  a  high  kinetic 
energy  in  some  regions  that  make  the  wave  function  oscillate  rapidly.  In  tunneling  regions,  the 
wave  functions  are  growing  or  decaying  exponentially  and  the  Schrodinger  equation  for  those 
regions  is  of  the  stiff  type.  To  obviate  this  difficulty,  we  divide  space  into  three  regions,  the 
left  side  for  which  energy  is  greater  than  the  potential,  middle  for  which  energy  is  less  than  the 
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potential  and  right  for  which  energy  is  again  greater.  (For  electrons  with  positive  kinetic 
energy  everywhere,  there  is  no  middle  region.)  With  two  linearly  independent  solutions  in  all 
three  regions,  a  solution  valid  everywhere  is  obtained  by  matching  the  solution  first  at  one 
turning  point,  then  at  the  other.  At  the  far  edge  of  the  left  and  right  sides,  the  potential  is 
assumed  constant  so  that  y(z)  =  exp  (±ikz).  For  the  propagating  regions,  the  wave  function  is 
recast  as  amplitude  and  phase, 

y(z)  =  R(z)  exp  [±i<j>(z)] 
and  the  Schrodinger  equation  becomes 

R"(z)  =  [v^z)  -  E]R  +  E/R3 

<J>'(z)  =  VE/R2  . 

In  the  tunneling  regions,  y  is  recast  reflect  its  exponentially  varying  behavior: 

y(z)  =  exp  [<()(/.)] 

and  the  Schrodinger  equation  becomes: 

<T(z)  =  ^)'(z)2  +  V(z)  -  E  . 

The  two  linearly  independent  solutions  can  be  obtained  by  starting  with  any  different  set  of  ini¬ 
tial  conditions.  In  practice,  the  equation  is  integrated  from  the  left  to  the  right  with  one  set  of 
initial  conditions  and  the  right  to  the  left  with  another,  this  avoids  the  problem  of  solving  stiff 
equations.  (It  is  also  possible  to  solve  the  Schrodinger  equation  in  a  WKB  approximation; 
however,  for  arbitrary  potentials  it  is  necessary  to  solve  the  equations  numerically  anyway,  and 
the  WKB  may  not  be  accurate  near  the  surface  where  the  potential  varies  rapidly:  it  is  there 
that  the  image  force  lowering  potential  the  MIGS  potential  (or  whatever  potential  is  responsible 
for  pinning  the  Schottky  barriers)  are  present.) 

It  is  well  known  that  electrons  traversing  the  barrier  can  tunnel  with  kinetic  energy  less 
than  the  barrier  height  (the  majority  of  current  is  tunneling  current  is  moderately  and  heavily 
doped  devices),  and  reflect  off  the  barrier  even  with  positive  kinetic  energy  [see,  for  example, 
Landau  and  Liftshitz  Section  50).  When  calculating  this  effect  in  Schottky  barriers,  it  has  been 
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customary  to  approximate  the  potential  barrier’s  shape  either  by  a  triangle  or,  for  reflection 
over  the  top  of  the  barrier,  by  a  parabola.  Our  program  permits  us  to  calculate  numerically  the 
true  reflection  coefficient;  it  is  found  to  differ  significantly  from  the  model  potentials.  Figure  1 
compares  the  reflection  for  a  potential  obtained  by  a  jellium  with  an  image  force  superimposed 
to  that  with  parabolic  potential  whose  curvature  matches  the  former  potential.  These 
differences  are  important  when  a  quantitative  comparison  with  experiment  is  desired.  In  cases 
when  the  band  structure  deviates  significantly  from  effective  mass  behavior  (this  is  especially 
important  for  tunneling  electrons),  one  can  easily  include  this  effect  using  a  two-band  model  or 
approximately  by  permitting  the  mass  to  be  position-dependent. 
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FIGURE  1  REFLECTION  PROBABILITY  OF  AN  ELECTRON  SCATTERED 
BY  A  SCHOTTKY  BARRIER 

In  summary,  development  of  the  tools  required  for  a  quantitative  study  is  proceeding  on 
schedule.  Within  the  next  few  months,  we  expect  to  be  able  to  assemble  the  pieces  to  make 
careful  studies  of  at  least  some  aspects  of  transport  through  a  Schottky  barriers.  Eventually  it 
is  expected  that  the  Schottky  barriers  work  also  will  evolve  into  an  effective  modeling  tool  for 
the  study  of  junctions  and  high-speed  electronics  devices. 
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The  next  section  presents  a  draft  of  a  paper  in  preparation  on  the  electrostatic  potential 
from  dopants  in  the  depletion  region.  As  mentioned  previously,  the  resulting  barrier  is  shown 
to  be  remarkably  similar  to  the  quadratic  potential  of  a  jellium,  as  is  usually  assumed.  The 
difference  between  these  potentials  acts  as  a  perturbation  to  electrons  traversing  the  barrier;  it 
is  shown  to  be  small  and  the  scattering  weak. 


5 


n  EFFECT  OF  DISCRETE  DOPANTS  IN  SCHOTTKY  BARRIERS 


Traditionally  the  barrier  in  Schottky  barriers  and  heterojunctions  is  modeled  by  a  one¬ 
dimensional  quadratic  potential,  as  would  obtain  from  jellium  of  charge  in  the  depletion  region. 
Fluctuations  in  the  potential  act  as  a  perturbation  to  electrons  traversing  the  barrier;  particularly 
important  are  the  fluctuations  in  the  plane  normal  to  the  barrier. 

At  first  glance  it  might  seem  that  the  difference  in  electrostatic  potential  from  a  random 
distribution  of  ionized  dopants  and  a  jellium  is  large  and  that  it  would  strongly  scatter  electrons 
traversing  the  depletion  region.  Boudvillc  and  McGill  (1985)  calculated  transport  using  a  sim¬ 
ple  model  for  the  coulomb  potential  and  found  the  effect  to  be  large.  This  section  presents  a 
more  complete  study  of  this  perturbation.  We  do  so  with  an  approximate  potential  and  later 
show  the  results  of  a  more  complete  Madelung  calculation  using  a  collection  of  100  randomly 
distributed  dopants. 

To  make  an  approximate  perturbation  to  the  true  potential,  consider  a  single  sphere  of 
radius  rs  and  volume  1/Nd,  where  Nd  is  the  dopant  density.  The  charge  density  inside  the 
sphere  is  made  out  of  a  point  charge  at  the  center,  compensated  by  uniform  background  of 
density  Nd.  The  potential  from  all  dopants  is  obtained  from  a  superposition  of  these  sphere 
potentials;  this  deviates  from  the  true  perturbation  because  the  charge  density  is  doubly  counted 
in  regions  of  overlapping  spheres  and  not  counted  in  void  regions.  The  net  charge  of  both  the 
true  and  approximate  perturbations  is  zero.  (There  is  an  additional  approximation  of  the  same 
order,  that  the  potential  from  the  superposition  of  charge  densities  is  taken  to  be  equal  to  the 
superposition  of  potentials  from  the  densities  separately — a  correction  comes  from  the  regions 
where  the  densities  overlap.)  This  approximation  essentially  the  same  as  the  atomic  spheres 
approximation  in  electronic  structure  calculations.  For  close  packed  lattices,  it  is  an  excellent 
one,  introducing  errors  of  few  percent.  For  a  random  distribution  of  dopants,  the  approxima¬ 
tion  should  still  be  reasonable.  The  potential  for  a  single  sphere  is,  in  atomic  Rydberg  units; 


(1) 


6 


for  r  <  rs  and  zero  otherwise.  (The  potential  must  also  be  reduced  by  the  dielectric  constant; 
this  reduction  will  be  implicit  in  all  expressions  for  the  potential  presented  in  this  paper.) 

TV  average  perturbation. 


♦  =  f'df ‘I'rfjlri  =  ^5- 


4fNd 


1/3 


(2) 


increases  as  Nd1/3,  much  as  the  image-force  correction,  varying  as  Nd/4.  This  approximate  <5>  is 
independent  of  the  distribution  of  dopants.  Its  magnitude  is  small  and  suggests  that  the  pertur¬ 
bation  is  weak.  For  example,  for  Nd  =  8  x  1018  cm-3  and  e  =  10,  <}>  =  15  meV,  barely  observ¬ 
able  even  by  the  best  IV  measurements  (appearing  a  correction  to  the  barrier  height)  and  small 
compared  to  the  image-force  lowering  potential. 

Now  consider  scattering  by  this  potential  when  dopants  arc  ordered  on  a  lattice.  For  any 
lattice  it  is  of  course  possible  to  obtain  the  potential  by  Ewald  summation,  but  our  approximate 
potential  should  be  a  quite  adequate.  When  there  is  a  lattice,  the  perturbation  is  most  informa¬ 
tively  expressed  as  a  Fourier  scries,  since  the  Fourier  components  are  the  “oscillator 
strengths,”  a  measuring  of  the  coupling  strength  of  one  state  to  another,  and  thus  the  scattering 
rate  between  states.  Components  of  the  approximate  perturbation  are  <J>(G)  with 

£<$,(?+ R)  =  NdX<D(G)  exp  iS  •  ?  .  (3) 

I i  cl 


Here  and  Cj  arc  real-  and  reciprocal-space  lattice  vectors,  respectively,  and  4>(q)  is  the 
Fourier  transform  of  4>(r), 

Nj<Kq)  =  ~  x  — ~T  Kqrs)3  +  3qrscos(qrs)  -  3sin(qrs)]  .  (4) 

-S|\  (qrs» 

Function  Nd$(q)  falls  smoothly  to  zero  from  its  maximum  value  of  3/5rs  at  q  =  0,  with  a  half- 
width  at  qrs  =  5.  Again,  the  strength  of  the  perturbation  is  small  (being  strongest  at  q  =  0), 
although  in  the  highly  doped  case  not  negligible  in  comparison  to  kT.  The  half-width  qrs  =  5 
is  a  measure  of  the  maximum  of  change  in  wave  vector  qf  an  electron  will  suffer  when  it 
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scatters.  For  example,  in  the  case  m  *  =  0.1  m,  Nd  =  8  x  1018  cm-3  and  e  =  10,  the  “oscilla¬ 
tor  strength”  is  15  meV  at  q  =  0  and  H2( 5/R)2/2m*  35  meV,  so  that  scattering  events  are  lim¬ 
ited  to  a  change  in  energy  of  approximately  35  meV. 

In  the  case  of  a  random  distribution  of  dopants,  the  same  gross  behavior  should  be  found 
as  before;  in  particular,  according  to  the  approximations  discussed  above,  the  average  perturba¬ 
tion  <j>  is  independent  of  dopant  distribution.  One  might  expect  that  the  random  distribution 
will  enhance  fluctuations,  because  two  or  more  dopants  may  cluster  together  and  combine  to 
form  a  local  potential  well  that  strongly  influences  scattering.  As  a  check  of  the  simple  model 
potential  described  above,  and  also  to  study  the  effect  of  randomness,  a  “sample”  depletion 
region  was  made  out  of  100  randomly  distributed  dopants.  The  dopants  were  placed  randomly 
in  the  bottom  half  of  a  supercell  of  width  w  and  length  2L;  the  top  half  contained  no  dopants 
and  thus  constituted  a  neutral  region.  Periodic  boundary  conditions  were  imposed  so  that  the 
suf>;rccll  repeated  itself  throughout  space.  The  electrostatic  potential  from  this  distribution  was 
obtained  by  Ewald  summation,  tabulated  on  a  mesh  of  50  x  50  x  50  points.  A  dopant  density 
Nd  =  (100  au)-3  =  8  x  1018  cm-3  and  a  length  L  =  332  au  was  chosen,  making  a  barrier  height 
of  1.39  Ryd  and  a  width  of  w  =  548  au.  (For  e  =  12,  this  corresponds  to  about  1.5  cV.) 

In  the  illustrations  shown,  the  z  axis  is  normal  to  the  depletion  region.  Figure  2  shows 
contour  plots  in  the  yz  plane  for  three  different  averages  of  the  potential  in  x.  In  Figure  2(a), 
the  contours  arc  shown  for  the  25th  plane  x,  in  Figure  2(b),  the  contours  are  averaged  in  x 
between  the  21st  and  30th  plane,  and  in  Figure  2(c),  the  contours  are  averaged  in  all  50  planes 
of  x.  As  expected,  the  potential  increases  approximately  parabolically  in  the  depletion  region 
(between  z  =  0  and  z  =  L),  and  the  fluctuations  die  down  as  more  planes  are  averaged.  Figure 
2(a)  is  not  very  meaningful  because  electrons  of  wave  vector  k,  will  be  smeared  out  in  x  with 
a  wavelength  27t/kx.  The  10  atomic  layers  of  Figure  2(b)  amount  to  about  50  A,  which 
corresponds  to  a  transverse  kinetic  energy  of  15  meV  for  m*/m  =  0.1.  Figure  2(c)  is  hardly 
distinguishable  from  a  parabolic  barrier  and  the  barrier  height  is  very  close  to  the  value  of  1.39 
Ryd  that  a  jellium  would  have. 

Figure  3  shows  the  deviation  from  a  jellium,  averaged  in  10  and  50  planes  in  x.  The  10- 
plane  average  is  seen  to  be  about  three  times  stronger  than  the  50-planc  average.  Figure  3(c) 
shows  the  deviation  in  potential  averaged  in  both  x  and  y.  It  is  seen  to  rise  smoothly  from 
approximately  zero  near  the  edge  of  the  depiction  region.  The  average  potential  in  the  deple¬ 
tion  region  can  visually  be  seen  to  be  close  to  10  mRy  the  approximate  potential  estimates. 
Figure  4  shows  the  Fourier  transform  of  the  potential  in  the  xy  plane,  for  z  in  at  the  edge  of 
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(al  CONTOURS  IN  THE  yz.  TAKEN  FOR  A  SINGLE  PLANE  IN  x 


0  2 
Of 

1  3 


lb)  CONTOURS  IN  THE  yz,  AVERAGED  IN  10  PLANES  IN  x 


(c)  CONTOURS  IN  THE  yz,  AVERAGED  IN  ALL  50  PLANES  IN  x 


FIGURE  2  CONTOUR  PLOT  OF  A  SCHOTTKY  BARRIER: 

CONTOURS  OF  CONSTANT  BARRIER  HEIGHT 
FOR  THREE  DIFFERENT  AVERAGES  OF  THE 
BARRIER  IN  x 
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(a)  CONTOURS  OF  CONSTANT  DIFFERENCE,  ANALOGOUS  TO  FIGURE  2<a| 


(c)  DIFFERENCE  WHEN  THE  POTENTIAL  IS  AVERAGED 
IN  THE  ENTIRE  PLANE  NORMAL  TO  THE  SCHOTTKY 
BARRIERS 


FIGURE  3  DIFFERENCE  IN  POTENTIAL  BETWEEN  THE  JELLIUM 
AND  THE  RANDOM  DISTRIBUTION  OF  DOPANTS 
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FIGURE  4  FOURIER  TRANSFORM  OF  THE  SCHOTTKY  BARRIER  POTENTIAL 
AT  THE  EDGE  OF  THE  DEPLETION  REGION 

Units  of  k  are  50/w,  where  w  =  548  au. 


the  depiction  region  and  in  the  center.  The  Fourier  transform  $(q,)  dies  out  very  rapidly  away 
from  q,  =  0,  indeed  too  rapidly  to  be  resolved  on  the  discrete  mesh  used.  This  rapid  decay  is 
perhaps  most  easily  interpreted  from  the  point  of  view  of  the  model  potential.  Taking  its 
Fourier  transform  over  a  random  distribution  of  dopants, 

Nd  fd3re^  '  =  Nd<Kq)  I  e1^ '  g  ,  (5) 

it  if 


the  sum  over  a  distribution  of  phases  ^  S  now  washes  out  all  of  the  Fourier  components  of 
this  approximate  perturbation  except  for  q  =  0. 

Thus,  it  appears  that  a  random  distribution  of  dopants  scatters  even  more  weakly  than  an 
ordered  lattice.  The  potential  in  either  case  appears  to  electrons  transversing  the  depletion 
region  to  be  remarkably  like  a  simple  quadratic  potential.  The  key  to  this  perhaps  surprising 
result  lies  in  the  long  range  of  the  coulomb  potential:  a  superposition  of  potentials  falling  off  as 
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1/r  combine  to  form  a  very  smooth  potential.  (The  singularity  in  a  1/r  behavior  at  short  range 
is  not  significant,  because  it  is  only  J(l/r)d3r  over  some  volume  characteristic  of  the  size  of  an 
electron  wavelength  that  matters.)  No  explicit  calculations  of  scattering  through  the  depletion 
region  were  carried  out  here,  although  this  is  rather  easy  to  do.  It  is  clear  from  the  discussion, 
however,  that  the  net  scattering  from  these  dopants  are  quite  small  (especially  in  the  more  phy¬ 
sically  significant  random  case),  particularly  when  averaged  over  a  distribution  of  electrons. 
Such  a  conclusion  is  the  opposite  of  the  conclusion  Boudville  and  McGill  (1985)  drew.  He 
did  carry  out  transport  calculations,  but  used  a  crude  potential  in  which  only  a  single  dopant 
was  present,  and  it  is  probably  in  the  differing  treatment  of  the  potential  that  our  conclusions 
differ. 
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m  AB  INITIO  ELECTRONIC  STRUCTURE  WORK 


Most  properties  of  interest  in  the  solid  are  governed  by  the  behavior  of  the  electrons,  and 
these  properties  can  (in  principle)  be  calculated  if  the  underlying  Schrodinger  equation  (or  the 
relativistic  extensions  to  it)  could  be  solved.  In  particular,  the  heat  of  formation  of  metal  ada¬ 
toms  on  the  surface  of  semiconductors,  the  interfacial  dipole  that  pins  most  Schottky  barriers, 
and  of  course  the  energy  bands  are  the  kinds  of  quantities  of  particular  relevance  to  Schottky 
barriers.  The  greatest  portion  of  our  efforts  in  the  past  year  have  been  devoted  towards 
development  of  an  efficient  method  for  making  such  electronic  structure  calculations,  one  that 
is  sufficiently  fast  that  it  can  be  applied  to  large  systems  (by  ab  initio  standards)  containing 
several  dozens  of  atoms. 

While  it  is  nearly  impossible  to  solve  that  equation  exactly  for  a  real  system,  the  local- 
density  approximation  to  it  has  been  found  to  yield  remarkably  good  results  under  a  wide 
variety  of  circumstances.  Hohcnbcrg  and  Kohn  proved  long  ago  that  the  Schrodinger  equation 
is  a  functional  only  of  the  electron  density;  the  functional,  however,  is  unknown.  The  local- 
density  approximation  is  a  simplification  of  the  unknown  density  functional,  and  it  is  possible 
(although  difficult)  to  make  systematic  corrections  to  it.  It  is  ab  initio,  in  the  sense  that  no 
adjustable  or  empirical  parameters  enter  into  the  theory.  Extensive  experience  in  a  truly 
remarkable  range  of  applications  has  generated  a  high  degree  of  confidence  in  the  LDA,  and  it 
is  widely  believed  that  the  it  is  sufficiently  accurate  to  predict  a  broad  range  of  properties,  par¬ 
ticularly  mechanical  properties.  It  can  be  equally  well  applied  to  any  element  in  the  periodic 
table  and  to  any  arbitrarily  complicated  system,  barring  formidable  difficulties  in  accurately 
obtaining  solutions  to  the  LDA  for  large  systems. 

Its  greatest  drawback  is  that,  at  least  as  traditionally  applied,  it  is  quite  complicated.  First 
a  basis  set  must  be  chosen  in  which  matrix  elements  of  the  Hamiltonian  can  be  calculated. 
Then  the  eigenvalue  problem  must  be  solved  to  generate  a  sum  of  one-electron  energies  and  a 
charge  density;  next,  Poisson’s  equation  must  be  solved  to  evaluate  the  Hartrcc  potential  and 
the  electrostatic  energy.  Finally,  the  total  energy  is  obtained  from  the  eigenvalue  sum,  plus  a 
double  counting  correction  in  the  Hartrce  potential.  Ab  initio  methods,  as  least  in  their  tradi¬ 
tional  application,  additionally  require  self-consistency  in  the  potential,  i.c.,  the  potential  as 
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calculated  from  the  eigenstates  of  the  Hamiltonian  must  be  identical  to  the  potential  that  gen¬ 
erated  the  Hamiltonian  in  the  first  place.  Unt:  recently  as  last  year,  the  available  computa¬ 

tional  techniques  were  still  regarded  as  far  too  cumbersome  to  calculate  properties  of  anything 
but  very  small  systems,  e.g.,  systems  with  20  atoms  or  less.  The  bulk  of  our  effort  this  year 
has  gone  into  development  of  an  efficient  method  that  solve  the  local  density  equations  far 
more  efficiently  than  previously.  This  opens  many  doors  that  were  previously  closed  to  ab  ini¬ 
tio  calculations;  in  particular,  it  should  be  possible  to  make  ab  initio  studies  of  the  metal- 
semiconductor  interface. 

An  accurate  solution  of  the  local  density-functional  is  a  complex  and  computationally 
intensive  task.  This  has  until  now  been  particularly  so  with  techniques  that  attempt  to  solve 
the  local-density  equations  with  essentially  no  approximations  (the  APW,  LAPW,  and  pseudo¬ 
potential  methods).  Calculations  using  smaller,  more  efficient  basis  sets,  in  particular  the  linear 
muffin  tin  orbitals  (LMTO)  method,  traditionally  make  simplifying  approximations  to  the 
potential.  The  LMTO  method  is  far  more  efficient  than  its  sister  methods,  but  the  approxima¬ 
tions  to  the  potential  render  it  suitable  only  in  conditions  of  high  symmetry. 

The  bulk  of  our  efforts  this  year  have  been  spent  in  collaboration  with  Dr.  Methfcssel  at 
the  Max  Planck  Institut  fur  Festkorperforschung  in  Stuttgart,  to  develop  a  new  LMTO  method 
that  removes  the  approximations  to  the  potential.  Because  the  method  is  only  now  being  com¬ 
pleted,  very  few  results  are  available.  However,  the  early  results  are  very  promising.  The 
LMTO  method,  without  any  shape  approximations  to  the  potential,  and  with  an  enlargement 
over  the  conventional  basis  can  be  as  accurate  as  the  computationally  intensive  methods  (Meth¬ 
fcssel,  1988;  Methfcssel  and  van  Schilfgaarde,  in  preparation).  Silicon  has  been  used  as  a  test¬ 
ing  ground;  it  makes  for  a  particularly  stringent  test  for  the  LMTO  method  because  of  its  open 
tetrahedral  structure  (the  LMTO  method  prefers  close-packed  structures).  Table  1  attests  to  the 
remarkable  precision  of  both  density-functional  theory  and  this  new  LMTO  method. 

This  basis  set  of  22  orbitals/atom  is  sufficient  to  solve  the  local-density  equations  to  an 
absolute  precision  of  about  1  mRy.  Comparable  accuracy  using  the  pseudopotential  method 
requires  about  1000  plane  waves  per  atom  (similarly  for  the  LAPW  method).  Because  the 
computation  time  for  a  given  calculation  of  the  energy  eigenvalues  increases  as  the  third  power 
of  the  number  of  orbitals  in  the  system,  it  is  obvious  that  this  method  is  orders  of  magnitude 
faster  than  the  pseudopotential  method  for  comparable  accuracy.  The  reason  why  the  LMTO 
method  is  so  efficient  is  that  its  orbitals  are  tailored  to  the  potential  of  the  system,  and  so  very 
rapid  convergence  is  obtained. 


14 


Table  1 


PROPERTIES  OF  SILICON  CALCULATED 
SELF-CONSISTENTLY  FROM  22-ORBITAL  BASE, 
FULL  POTENTIAL 

Data  from  Methfessel  (in  preparation) 


Parameter 

Theory 

Lattice  constant  (a.u.) 

10.26 

10.23 

Cohesive  energy  (eV)* 

4.8 

5.23 

Bulk  modulus  (Mbar) 

0.98 

0.987 

cu  -  c12  (Mbar) 

1.02 

1.03 

c44  (Mbar) 

0.80 

0.83 

Phonons: 

TO(T)  (THz) 

15.53 

15.52 

Kxyz  (eV/A) 

-35.1 

39.1 

TO(X)  (THz) 

13.90 

13.75 

LAO(X)  (THz) 

12.32 

11.82 

TA(X)  (THz) 

4.49 

4.50 

Griineisen  parameters: 

TO(T) 

0.98 

TO(X) 

1.5 

1.51 

LAO(X) 

0.9T 

1.03 

TA(X) 

1.4 

1.42 

This  is  essentially  the  result  of  local-density  theory,  and 
the  error  is  almost  entire  owing  to  failure  of  the  local 
density  in  the  free  atom,  as  opposed  to  the  solid. 

IThere  is  a  large  uncertainty  in  the  measurement  of  this 
quantity  (Manuel  Cardona,  private  communication). 


LMTO  calculations  have  been  performed  on  large  systems,  including  the  NiSi2-Si  inter¬ 
face,  using  nine  interfacial  layers  comprising  80  atoms  (Das  et  al„  submitted;  Das  et  al.,  in 
press)  and  one  on  a  216-atom  silicon  cluster  thought  to  resemble  amorphous  silicon.  This 
demonstrates  the  feasibility  of  applying  the  LMTO  method  to  large  systems;  our  present  task  is 
to  do  the  same  with  the  accurate,  full-potential  version  of  the  LMTO  method. 
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Another  theoretical  development  important  to  ab  initio  treatment  of  large  systems  comes 
in  new  approaches  to  the  density-functional.  Self-consistency  can  be  a  real  obstacle  in  conven¬ 
tional  applications,  first  because  it  is  slow,  and  secondly  because  it  is  hard  to  reach,  owing  to 
the  many  more  degrees  of  freedom.  A  recent  development  by  Car  and  Parrinello  unites 
density-functional  theory  and  molecular  dynamics,  and  self-consistency  in  the  electronic  struc¬ 
ture  is  achieved  simultaneously  with  the  motion  of  the  nuclei  (Car  and  Parrinello,  1985).  This 
has  been  successfully  applied  to  a  (001)  twist  boundary  in  germanium,  using  the  pseudopoten¬ 
tial  method  a  local  approximation  to  the  exact  LDA  pseudopotential  (Payne,  Bristowe,  and 
Joannopoulos,  1985). 

Another  recently  developed  technique  (Harris,  1985;  Foulkes,  1987)  is  even  simpler  than 
the  above  mentioned  Car  and  Parrinello  method,  and  more  efficient  for  the  LMTO  method. 
These  techniques  are  derivatives  of  a  method  originally  published  by  Harris  (1985),  and 
independently  and  more  completely  by  Foulkes  (1987).  The  essential  point  here  is  to  exploit 
the  variational  property  of  the  Hamiltonian,  because  of  which  errors  in  the  total  energy  are 
second  order  in  the  difference  between  the  guessed  input  potential  and  the  self-consistent 
potential.  Rather  than  carrying  calculations  to  self-consistency,  one  attempts  to  construct  an 
input  potential  sufficiently  close  the  the  self-consistent  one  as  to  render  unnecessary  any  steps 
to  self-consistency.  Another  key  in  this  technique  is  that  the  density-functional  can  be  written 
in  another  form  so  as  to  require  only  a  guessed  input  potential  and  the  output  band  structure 
energy. 

The  Hohenbcrg-Kohn  density  functional  at  a  density  n  is 

E[n]  =  T[n]  +  E„artree[n]  +  E*c[n]  +  JnVen[n]  , 

where  EIIartree[n)  is  the  electrostatic  energy  l/2j  n(r)  n(ri)/ 1  r  -  rild3rdV  of  the  electrons,  jnVc 
is  the  energy  of  the  electrons  interacting  with  the  nuclei  and  Exc[n]  is  the  exchange-correlation 
energy. 

The  kinetic  energy  T  is  not  directly  calculable,  but  it  can  be  obtained  by  solving  the 
Schrbingcr  equation  for  its  eigenstates  y,  and  eigenvalues 

<T  +  vm)  ^  ^  . 
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where  T  =  -  1/2V  2  and  vm  is  a  guessed  input  potential.  Solution  of  the  one-electron  equa¬ 
tions  for  v| f{  yields  an  output  charge  density 

"out  =  Si  v‘  Vi  . 

where  is  the  sum  over  occupied  orbitals.  Were  the  self-consistent  potential,  n^,  would 
be  the  self-consistent  density  no  as  well,  and  vm  would  then  be 

vm  =  Ven(r)[no]  +  J  -^2-  dV  +  mc[no]  , 

but  (and  this  is  the  central  point  of  this  approach)  need  not  be  calculated  from  any  density 
(Foulkes,  1984).*  Because  of  the  variational  principle,  any  guessed  input  potential  incurs  errors 
of  order  (nou,  -  no)2.  Solution  of  the  above  equations  yields  an  expression  for  the  kinetic 
energy: 

TfaoJ  =  £  J  V*  TVi  =  E  £i  -  \  Vin"out  • 
i  i 

Making  a  functional  Taylor  series  in  T, 

T[nJ  =  T[noul]  +  J  (n,n  -  noul)(5T/8nin)  +  0(nm  -  i^J2  , 

and  using  8T/5nin  =  -vlrl  +  constant,  Foulkes  (1987)  obtained  a  new  expression  for  the  total 
energy 

£[ttin>  "h  J  (^en— 

i 

4"  f^HartredTin]  E^fn^]  +  0(11^1  —  n^)  +  0(noul  —  n@) 

This  last  expression  differs  from  the  self-consistent  one  by  errors  of  second  order  in  the  den¬ 
sity.  It  is  exactly  equivalent  to  the  density-functional  as  originally  formulated,  but  is  amenable 
to  approximation,  especially  with  respect  to  the  LMTO  method.  In  particular,  the  LMTO 


*The  idea  that  the  input  potential  v  need  not  be  derived  from  an  input  density,  also  due  to  Foulkes,  is 
recent  and  is  as  yet  unpublished. 
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method  considerably  simplifies  if  it  is  sufficient  to  construct  a  starting  potential  from  a  super¬ 
position  of  spherical  atom-centered  potentials.  M.  Methfessel  (in  preparation)  has  constructed 
a  test  potential  in  silicon  and  found  good  agreement  with  the  fully  self-consistent  results. 
Importantly,  the  same  test  potential  can  be  put  in  different  structures  (e.g.,  beta-tin,  diamond, 
simple  cubic),  so  it  is  transferable  to  new  environments.  Self-consistency  is  thus  obviated. 
This  method  is  in  practice  essentially  as  efficient  than  the  semiempirical  tight-binding  method, 
but  has  the  full  support  and  precision  of  density-functional  theory.  Indeed  it  very  strongly 
resembles  the  tight-binding  method  since  essentially  all  that  enters  are  the  sum  of  one-electron 
energies.  One  of  the  greatest  achievements  of  the  Harris  functional  is  to  make  a  formal 
justification  of  most  of  the  ansatz  used  implicitly  in  semiempirical  tight-binding  methods.  The 
method  outlined  here  is  far  more  efficient  than  the  Car-Parrinello  method  and  is  well  suited  for 
large  systems. 

A  related  approach  is  to  make  a  guessed  potential  out  of  a  superposition  of  free-atomic 
charge  densities.  This  is  more  complicated  to  implement  since  the  potential  is  no  longer  a 
superposition  of  spherical  potentials,  but  the  results  seem  about  as  good  (Polatoglou  and  Meth¬ 
fessel,  1988)  except  for  very  ionic  systems  such  as  NaCl. 
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IV  SUMMARY  AND  FUTURE  WORK 


One  step  in  the  current  suite  of  LMTO  programs  involves  a  technique  that  solves 
Poisson’s  equation  and  the  exchange-correlation  potential  and  energy  in  the  interstitial  (Meth- 
fessel,  1988;  Methfessel  and  van  Schilfgaarde,  in  preparation).  The  idea  is  to  fit  the  density  as 
evaluated  at  the  surfaces  of  the  spheres  to  a  linear  combination  of  spherical  Hankel  functions; 
since  Hankel  functions  are  eigenfunctions  of  the  Poisson  equation,  the  Hartree  potential  is  trivi¬ 
ally  obtained  as  a  linear  combination  of  Hankel  functions  once  the  fit  to  the  density  is  known. 
While  ingenious,  the  method  as  originally  designed  was  not  suitable  for  large  systems  because 
of  the  method  used  to  fit  the  density  itself.  However,  we  have  recently  shown  (Methfessel  and 
van  Schilfgaarde,  in  preparation)  that  procedure  is  mathematically  equivalent  to  a  generalization 
of  Andersen’s  tight-binding  transformation  (Andersen,  1985),  and  that  the  procedure  can  be 
made  far  more  efficient  than  previously.  In  particular,  it  is  possible  to  solve  Poisson’s  equation 
in  real  space;  the  method  is  equally  well  suited  to  molecules.  It  was  principally  in  this  step 
that  the  method  as  originally  implemented  was  not  suitable  for  large  systems;  this  difficulty  is 
being  removed  with  the  current  work. 

A  number  of  improvements  need  be  made  to  make  it  fully  operational,  for  example,  the 
full-potential  version  of  the  code  is  as  yet  neither  relativistic  nor  spin  polarized.  Without  the 
atomic  spheres  approximation  (in  which  the  potential  is  constructed  out  of  large  overlapping 
spheres),  semicore  states  are  a  greater  problem  than  previously  and  probably  a  two-panel  facil¬ 
ity  will  be  required.  These  are  a  number  of  details  that  need  to  be  filled  out  for  the  method  to 
be  generally  applicable  to  any  system.  We  believe  that  this  new  method  shows  a  great  deal  of 
promise  and  will  find  wide  application  in  the  years  to  come.  Because  it  is  both  very  accurate 
and  very  fast  we  believe  it  will  ultimately  displace  all  other  present-day  ab  initio  techniques  for 
calculating  electronic  structure  as  the  method  of  choice.  For  example,  it  should  be  possible  to 
study  the  metal-semiconductor  interface,  and  the  role  of  impurities  there  and  interface  states 
there.  With  a  Cray  II,  it  is  possible  for  the  method  to  perform  calculations  on  several  hundreds 
of  atoms. 
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